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1. Introduction 

Compute power increases every year and allows us to reduce statistical errors of lattice sim¬ 
ulations. Systematic errors however cannot be cured as simple as that and may pose challenging 
problems. A well-known unknown is the uncertainty due to the finite lattice spacing, which may 
have a significant effect. For determinations of renormalization factors in the context of RI-MOM 
schemes [ 1 ] the finite lattice spacing causes for instance so called hypercubic lattice artifacts in 
the data for the renormalization factors. Their origin lies in the breaking of the Euclidean con¬ 
tinuum symmetry 0(4) to the hypercubic subgroup of reflections and permutations, H( 4). This 
has the consequence that lattice results for the renormalization factors in momentum space depend 
on the direction of momentum, i.e. the orbit under the action of H( 4). To leading order in the 
momentum this discretization effect typically grows, but their exact functional form depends on 
the operator, the lattice spacing and the momentum. For the renormalization scale /./ one therefore 
requires A 2 <C p 2 <C X ■ A is the typical scale where nonperturbative effects are present and so the 
first condition ensures that perturbation theory can be applied, e.g., to connect renormalization fac¬ 
tors to the MS scheme. The second condition ensures discretization effects arc small. In practice, 
however, both conditions are not met simultaneously and thus hypercubic artifacts often spoil the 
accuracy of the determination. 

There arc few ways to alleviate this problem. One is to consider only diagonal momenta 
ap = a(q.q.q.q), where q = 2 nk^/N^, k u = — N^/2+ 1,... ,N^/2 and /V u is the number of lattice 
sites in direction /./ = 1.2.3,4. In these momentum directions hypercubic artifacts are much smaller 
than in any other direction but as we will see a significant amount still remains. 

A more sophisticated way to remove these discretization errors is perturbative subtraction. 
There, the idea is (1) to calculate the hypercubic artifacts Sf/(a.p) for a lattice quantity f/j at (a. p) 
within lattice perturbation theory (LPT) and then (2) to subtract this from the nonperturbative data 
for (fi at . If this subtraction is complete the corrected ()) al depends on a 2 p 2 rather than ap: 

0corr(a 2 p 2 ) = @iat(a,p) ~ S0(a,p). ( 1 . 1 ) 

In practice, such a subtraction has been used in [2] using diagrammatic LPT calculations at 1- 
loop order in §5 and also in [3, 4] to order a 2 . However, a 1-loop subtraction often appears to be 
insufficient. 

We therefore explore Numerical Stochastic Perturbation Theory (NSPT) as a tool to calculate 
higher-loop corrections to 8(f(a. p). We have successfully pursued this approach already for the 
gluon and ghost propagators in Landau gauge [5]. Here we will show that NSPT works equally 
well for RI'-MOM renormalization factors. To make the hypercubic artifacts even stronger for our 
study we use standard (unimproved) Wilson fermions. It is also sufficient to look at quenched QCD 
to see that our method works. This choice was simply for practical reasons and unquenched sim¬ 
ulations do not impose any problems for future determinations. Neither does the use of improved 
actions: If one is able to subtract large discretization effects which go along with unimproved Wil¬ 
son fermions, one should have good chances to remove the smaller errors one typically sees for 
improved actions and operators. 
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2. RI'-MOM factors from Numerical Stochastic Perturbation Theory 


Numerical Stochastic Perturbation Theory (NSPT) makes use of the Langevin process to pro¬ 
vide expectation values for observables in the perturbative regime of lattice QCD. It is reviewed 
in detail in [6], and it suffices to mention here that for the Langevin updates and all measure¬ 
ments we treat link fields (NSPT configurations) and derived quantities as power expansions in 

(l- l/2 =g 0 /V2N c , 

U x .n = 1 + £ J3 - i/2 uil + 0($ ~"/ 2+1 ). (2.1) 

i= 1 

Hence, all algebraic operations are performed with respect to this expansion and so additions and 
multiplications are as for polynomials: [A + B]^ = A in> + B [n] and \AB] ( " : = re¬ 

spectively. An expansion can be inverted using the recursive formula 


[A-i] W = -[aW]" 1 (2.2) 

(=0 

The latter we use for example to calculate two- and threepoint functions from the standard Wilson 
Dirac operator, making use of the well-known zeroth-order inverse, the tree-level Feynman propa¬ 
gator So- In our setup, we solve the Langevin equation using the simplest Euler integration scheme 
for three different step sizes £ = 0.01,0.02,0.03 and extrapolate our results afterwards linearly to 
£ = 0. We simulate lattices of sizes N 4 = 16 4 ,24 4 and 32 4 to resolve finite size effects. 

Renormalization factors in the RI'-MOM scheme are defined in the chiral limit, for example, 
via the //(4)-invariant renormalization condition [3] 


7 ri'-mom~ 
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(2.3) 


where Z RI " MOM (/r) is the renormalization factor for a flavor non-singlet current in the RI'-MOM- 
scheme and Z RI _MOM is the wave function renormalization 

z Ri'-MOM fu) = 

v 12£ v sin 2 (ap v ) 

On the lattice, one defines the amputated vertex function, r(p) = S l (p)G(p)S l (p) using the 
threepoint function G(p) we get from the expectation value of a local quark current (Ax) = 
i7(x)r<i(.r), where I is an interpolator with the desired quantum numbers: 

n 12 

G(p) = — Y,exp(-ip-(x-y))(u(x)u(z)rd(z)d(y)). (2.5) 

' rv7 


(2.4) 


These expressions arc suitable both for the perturbative and the nonperturbative case and both we 
need here. So we proceed in a similar way: As fermionic n-point functions depend on the gauge, we 
first fix both the NSPT and the nonperturbative (quenched lattice QCD) configurations to Landau 
gauge and measure the two- and threepoint functions. The only difference is in the use of the 
algebraic operations for NSPT and how the chiral limit is achieved: The nonperturbative Z-factors 
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are obtained on a 32 4 lattice from a linear extrapolation to zero quark mass of data for three values 
of the hopping parameter K = 0.1489,0.1507 and 0.1520 [7]. For NSPT we can measure directly 
in the chiral limit, as we only need the tree-level Feynman propagator for the inversion of the Dirac 
operator (cf. (2.2)) to obtain all n-point functions. 


3. Hypercubic Corrections 


Our perturbative and nonperturbative lattice configurations were thermalized with respect to 
the same (quenched) action. Data accumulated on either set should therefore show similar hyper¬ 
cubic artifacts. For the perturbative data we also know - up to a constant - their exact values in the 
continuum limit. This allows us to construct the hypercubic artifacts to a given loop order n 


8Z , 


(») 


(ap)=Z^(ap)-z£ > act (ap), where Z, 


'RI'-MOM V 


A n ) 


t(ap ) = <7„,o + ^d„,/log((np)" 
;=1 


(3.1) 


and the r/, M ->o can be calculated from the anomalous dimension of the operator under consideration 1 
so that only the constant d n o remains to be fixed from the NSPT data. 

In practice, however, the data is still afflicted with finite volume effects and possibly large 
autocorrelation times from the Langevin process. Thus, we do fits to the form 
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where the first term is the finite constant we are ultimately interested in and the second term re¬ 
moves the logarithmic dependence. The third term parameterizes the finite size effects in a similar 
way as in [9] and the last sum is a suitable form to describe the hypercubic artifacts in our fit range 
[( a P)mi n - ( a P)max\ U P to a degree d, corresponding to a fj(a 2d ) removal. The free parameters of the 
fit are d n .o, cy and c,-j and their number and values depend on the choice of the fit range and d. The 
latter needs to be increased with the fit range to describe the finite lattice spacing effects. 

To estimate the error correctly, and to cover the interesting region of small ( ap) 2 while avoid¬ 
ing too strong systematic effects, we bootstrap over our data points and do various fits with different 
values of {ap) 2 nin , (ap)j nax and d. Finally, we take the weighted mean of the individual results using 
the probability density from the corresponding y 2 distribution. 

Our results for d n o can be read off from Table 1 a) up to three loop order. Table 1 b) com¬ 
pares the 1-loop results with values from the literature. Overall, the values are in good agreement, 
although the given statistical errors still seem to underestimate the true errors. 

With our values for d„fl, we can now quantify the hypercubic discretization errors, using 
Eqs. (3.1). The left panels of Fig. 1 show these corrections summed up for the vector and scalar 
renormalization factor setting /3 = 6.2. These corrections can then be used to improve the nonper¬ 
turbative data thermalized for the same /3, 


3 

Z corr (a,p) =Z NPR (a,p) - £/3~'<SZ (,) (a,p). 

i= 1 


(3.3) 


Data for Zy and Z$ before and after this subtraction is shown in the right panels of Fig. 1. 
'See, e.g., [8] for a similar calculation. 
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(a) Results for finite contributions (b) Comparison to LPT 


i 

4 , 

dfd 

4o 

4o 

NSPT 

LPT 

ref. 

i 

-1.045(3) 

-0.808(3) 

-0.862(3) 

-0.868(3) 


-1.045(3) 

-1.044510 

[ 10 ] 

2 

-1.93(2) 

-1.36(2) 

-1.65(3) 

-1.47(2) 

^AV 

-0.808(2) 

-0.800249 

[ 10 ] 

3 

-5.82(5) 

-3.79(3) 

-5.42(14) 

-4.06(4) 

7 ( 1 ) 

-0.862(3) 

-0.858819 

[ 11 ] 






^( 1 ) 

-0.868(3) 

-0.862146 

[ 10 ] 


Table 1: a) Final results for the finite contributions - i.e. the dip in Eq. (3.1) - of the renormalization 
constants for the vector (V), axialvector (AV), scalar (S) and tensor (T) currents, b) Comparison between the 
NSPT results of this work and known 1-loop LPT calculations. 


4. Boosted Perturbation Theory 


Tadpole contributions to perturbative expansion coefficients are often large in lattice perturba¬ 
tion theory. Such diagrams do not have a continuum analogue and are the main reason for the slow 
convergence of LPT [12]. Of course, also the NSPT results obtained here and the subtraction of 
discretization errors suffer from this problem. 

The solution is to rescale the expansion parameter by the plaquette [12]. The rationale for that 

is that tadpole diagrams always contain a purely gluonic loop which is proportional to the plaquette. 

2 

Hence, one defines a “boosted” or “tadpole-improved” coupling g^ —> g^ = or equivalently, 
/3 — > fa = Pj3, where P = 1 + £ i=1 [3 ‘p il> is measured perturbatively but in the usual way. Both 
the unboosted and boosted series of a quantity 0 with coefficients d^ and df, respectively, should 
converge to the same value: 

oo oo 

(4.1) 

i =0 1=0 

By inserting [3 1 = P [3 h 1 , we express the unimproved coupling in terms of the new, shifting the 
dependence on the old expansion parameter one order higher. To 3-loop order we find by comparing 
the coefficients of [3 b 1 


( 0 ) 

Ou = 


Oq\ 4 1} = °o’) °b ] = °o 2) +P (1) °o ) > °i 3) = °o } + 2p (l) o {2) +P ( 2 ) o ( 0 1) + (p (1) ) 2 o ( 0 l \ 


(4.2) 

so that one obtains the improved n-loop estimate 0 ~ L?=oftr*°i ’ which is, for any finite trunca¬ 
tion at order n, different from the original expansion. 

In our case, (3 = 6.2 gives P = 0.623, so that fa = 3.864. It turns out that boosted perturbation 
theory is very efficient in all cases: The boosted 1-loop subtraction has already a larger effect than 
the unboosted 3-loop subtraction. The right-hand sides of Fig. 1 show the subtractions to different 
loop-orders and compare the boosted with the unboosted results for the examples of the scalar and 
the vector renormalization factor. 

Boosted Perturbation Theory (BPT) has already been applied in the subtraction of 1-loop dis¬ 
cretization effects [ 2 ] and, as it turns out now, this is a very successful and efficient way to remove 
most of the artifacts. Nevertheless, also in boosted perturbation theory, the 3-loop subtraction gives 
in most cases a significant improvement over the lower orders. A pretty stable plateau is reached 
over nearly the whole range of momenta. The scalar renormalization factor seems to be an ex¬ 
ception to the rule and raises slight doubts on the applicability of boosted perturbation theory in 
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fop) 2 



fop) 2 


Figure 1: (Unboosted) Hypercubic corrections (left) and resulting plateaus for the unimproved (01), un¬ 
boosted improved (11 to 31) and boosted (11-bpt to 31-bpt) RGI vector (top) and scalar (bottom) renormaliza¬ 
tion factor. 
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Figure 2: Fits of the Vector (left) and Scalar (right) RGI renormalization factor to a constant. The fit range 
was fixed to a 2 p 2 £ [2,12]. The abscissa labels the degree of data improvement: 01 = unimproved, nl(-bpt) = 
improved by our corrections up to n-loop order in the bare (boosted) coupling. Errors are purely statistical. 


certain cases: The boosted 1-loop curve lies on top of the boosted 3-loop curve, whereas the 2-loop 
subtracted result lies even above the unboosted 1-loop subtraction. This can also be seen from 
Fig. 2 which shows fits to a smaller plateau region for the two examples, using different methods of 
improvement. The reason lies in the smallness of the (unboosted) 2- and 3-loop coefficients of the 
correction, so that the boosted 2- and 3-loop coefficients are dominated by the (negative) plaquette 
contributions (cf. Eq. (4.2)). This has the consequence that the boosted coefficients alternate in sign 
and the boosted subtraction probably alternates around the true value. Thus any finite subtraction 
might well “overshoot” the true value. 
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5. Conclusions 

Using NSPT we were able to calculate the leading hypercubic lattice corrections of several 
hadronic operators in quenched (boosted) lattice perturbation theory (LPT) with Nf = 2 valence 
Wilson fermions. Our results serve up to 3-loop order in an expansion of the gauge coupling 
gQ (resp. gjf and allow us to remove almost all hypercubic lattice artifacts in the corresponding 
nonperturbative data for the renormalization factors Zy, Zav, Zs and Zj in the RGI scheme. 

Looking at these corrections in more detail we find that a calculation at 3-loop order in g^ 
provides about the same corrections as a 1-loop calculation in gj r Moreover, in contrast to g$, the 
additional corrections from a 2- or 3-loop calculation in g^ would provide only a negligible further 
improvement for Zy, Zav, Zv and Zj (see Fig. 2), if diagonal lattice momenta are considered. 
This might be different for operators with derivatives and for off-diagonal momenta, which will 
be further analyzed in a forthcoming publication. Note also that boosted LPT leads to alternating 
contributions for Z<, for different loop orders. 

To provide a proof of principle for our approach, we chose a partially quenched setup with 
unimproved Wilson fermions, but we are not restricted to that. We have demonstrated that with 
NSPT one can estimate the momentum dependence for the leading hypercubic lattice corrections 
of different hadronic operators beyond the 1-loop order. Such calculations are indeed challenging 
but not as they were within traditional LPT, in particular with respect to all the improved actions 
used in state-of-the-art lattice QCD simulations. 
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